Poisson integrators for Volterra lattice equations 



T. Ergeng, B. Karasozen* 
Department of Mathematics and Institute of Applied Mathematics 
Middle East Technical University, 06531 Ankara- Turkey 



Abstract 

The Volterra lattice equations are completely integrable and possess bi-Hamiltonian struc- 
ture. They are integrated using partitioned Lobatto IIIA-B methods which preserve the 
Poisson structure. Modified equations are derived for the symplectic Euler and second or- 
der Lobatto IIIA-B method. Numerical results confirm preservation of the corresponding 
Hamiltonians, Casimirs, quadratic and cubic integrals in the long-term with different orders 
of accuracy. 
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1 Introduction 

The preservation of qualitative properties of discretized models of continuous systems became 
more and more important in recent years. The development of symplectic integrators for Hamil- 
tonian systems, construction of geometric integrators which preserve symmetries, reversing sym- 
metries and phase space volume of the underlying differential equations are some examples. For 
Hamiltonian systems in non-canonical form with a non-linear structure matrix, i.e Poisson sys- 
tems, there does not exist any general structure preserving integrator similar to the symplectic 
methods for canonical Hamiltonian systems (Sec. VIL2 ^Sj) . There are some well developed 
methods based on generating functions and Hamiltonian splitting for the Lie- Poisson systems, 
i.e. non-canonical Hamiltonian systems with a linear structure matrix (Sec. Vn.2.6 [S]). Pois- 
son systems arise especially as Hamiltonian pde's like Korteweg de Vries (KdV) and nonlinear 
^ Schrodinger equation with infinitely many integrals. Semi-discretization of these pde's in space 
by preserving the integrals and Hamiltonian structure results in integrable lattice equations; 
ode's with a Poisson structure and a finite number of integrals. 

After introducing the Volterra lattice equations in the next section, we discuss in Section 3 
symplectic integrators like the implicit mid-point and the symplectic Euler method applied to 
the well-known two dimensional Lotka- Volterra equation. Preservation of the Poisson structure 
of the Volterra lattice equation by symplectic Euler and Lobatto IIIA-B methods is proved and 
the corresponding modified equations are derived. The numerical results show the preservation 
of the conserved quantities by these methods. 



* e-mail: bulent@metu.edu.tr 



1 



2 Volterra lattice equation 



The m-dimensional Volterra lattice 

Vi = ViiVi+i - Vi-i), 



1, . . . ,m 



(1) 



for even m and with periodic boundary conditions i/m+i = Vi, i = 0,1, . . . and with yi > was 
studied first in 10 as an integrable system. It was shown that the Volterra lattice equation 
represents an integrable discretization of the KdV equation and of inviscid Burger's 

equation Besides these, the Volterra lattice equation describes many phenomena such as 
the vibrations of particles on lattices (Liouville model on the lattice) , waves in plasmas and the 
evolution of populations in a hierarchical system of competing species [Tl 01 ITU]. 
The system (1) possesses a bi-Hamiltonian structure [H I19j 



with respect to the quadratic and the cubic Poisson brackets 

{yi,yi+i}o = ViVi+i, 
{yi,yi+i}i = yiVi+iivi + Vi+i), {yi,yi+2}i = yiyi+iyi+2- 

The corresponding Hamiltonians are 

m ^ m 

Hi=^yi, Hq = -Y^ log(yi). 

i=l i=l 

and the structure matrices Jo(2/) corresponding to the quadratic Poisson bracket has the form 

/ yiy2 -yWrn \ 

: ■•• ■■• ■•• 



(2) 

(3) 
(4) 

(5) 



My) 



\ yiy-n 



-yi-Wi yiyi+i 



ym—iym 



/ 



If the Poisson brackets are compatible, i.e. the sum of {•,-}o + {t}i is again a Poisson 
bracket, then the bi-Hamiltonian system (1) has a finite number of functionally independent 
first integrals = 1, . . . ,m in involution i.e. {Ii,Ij}k = for A; = 0, 1 and i ^ j, i.e.. the 
Volterra lattice is completely integrable with respect to both brackets. 

The Hamiltonian Hq is a Casimir with respect to the Poisson bracket {•, - jo, ^-^-iHo, F}q = 
for any function F{y). 

The Volterra lattice (1) represents an integrable discretization of KdV equation jl9j 



du du d^u 



(6) 



which also possesses a bi-Hamiltonian structure and has infinitely many integrals. The first 
three integrals of © are [01 EI 



" ( 1 (du\^ 1 



di. 
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The corresponding conserved quantities of the Volterra lattice (1) are the quadratic and cubic 
integrals |7] 

iq = "^2^^ + yiyi+i^ h = ^^vf + yiyi+i{yi + yi+i + yi+2)- (7) 

i=l i=l 

In discretized form, the mass conservation Zi corresponds to the Hamiltonian Hi, the momentum 
and energy integrals T2, X3 of the KdV equation (5) correspond to the first integrals Iq and Ic 
of the Volterra lattice (1) respectively. 

The Volterra lattice is time-reversible and is closely connected to the Toda lattice, which is 
also completely integrable and has a tri-Hamiltonian structure [2], JH]- The transformation of 
variables ai = y2iy2i-i, h = y2i-i - y2i-2 in (1) gives the Toda lattice 

dj = aj(6j+i - bi = ai-ai-i, i = l,...,m/2 

with periodic boundary conditions 

ao = 0^/2 > bm/2+1 = h- 



3 Poisson integrators 

Geometric integrators for the Poisson systems 

y = J{y)VH (8) 

with a skew-symmetric non-constant structure matrix J{y) were studied recently in several 
papers. For recent surveys on Poisson integrators see (Sec. VII. 2, 8 and U). For related 
material on Poisson systems see (Sec. VII. 2, jH]) , (Ch. 10, 13J and (Ch. 6 & 7 .[TH]). 
The Poisson bracket {F, G} for two smooth functions F{y) and G{y) is defined by 

{F,G}iy) = VF{yfj{y)VG{y) (9) 

which is bilinear, skew-symmetric {G, F} = —{F, G} and satisfies the Leibniz' rule {{E-F, G} = 
E-{F,G} + F- {E, G}) as well as the Jacobi identity {{E, {F, G}} + {F, {G, E}} + {G, {E, F}}). 

The structure matrix J{y) of Poisson systems does not need to be invertible as for canonical 
Hamiltonian systems with J~^. All odd-dimensional skew-symmetric structure matrices J{y) 
are singular. The structure matrices Jo{y) and Ji{y) of the periodic Volterra lattice (1) are also 
singular. Such systems are called degenerate Poisson systems. 

The functions Cj(y) satisfying {Gi,H} = are called Casimirs or distinguished functions, 
which are first integrals whatever H(y) is. 

There are two characteristics for the flow ipt{y) of the Poisson system (8): 

• the flow (pt{y) of the differential equation (8) is a Poisson map, i.e. 

'ft{y)J{y)'ft{yV = J{ft{y)) 

where '^t{y) denotes the Jacobian of 99, 

• and it respects the Casimirs of J(?/), i.e. Gi{ipt{y)) = Const. 
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A numerical method i/n+i = 4'h{yn) is called a Poisson integrator for the structure matrix 
J{y) if the transformation y„ Un+i is a Poisson map that respects the Casimirs. The Casimirs 
should be preserved by the Poisson integrator. But the Casimir functions can be arbitrary, in 
case of the Volterra lattice they are in logarithmic form, therefore their conservation depends 
on the special structure of the problem. 

Because each Poisson system is distinguished by the structure matrix J{y), a method will be 
Poisson integrator only for a specific class of structure matrices. Therefore symplectic methods 
used for Hamiltonian systems can not be directly applied to Poisson systems. But some of 
the symplectic integrators can preserve certain Poisson structures. An example of this kind is 
the Poisson system resulting from the Ablowitz-Ladik integrable discretization of the nonlinear 
Schrodinger equation JHI which is preserved by the symplectic Euler method. 

The two-dimensional Lotka- Volterra equation was studied by several authors. The symplectic 
Euler and Stormer-Verlet methods preserve the Poisson structure of (10) whereas the implicit 
mid-point rule does not (see pp. 238 0). There are also some non-standard methods([T^. pTj ) 
which preserve the Poisson structure of the two-dimensional Lotka- Volterra equation. 

In the following we will apply the symplectic Euler method and Lobatto IIIA-B methods 
to the splitted form of the Volterra lattice equation (1). Both belong to the class of splitting 
methods which have been successfully used as geometric integrators in recent years. For a survey 
of the splitting methods see |Hj and . 

For the application of partitioned Runge-Kutta methods of Lobatto type, we split the equa- 
tion (1) into two parts 

Ui = Ui{vi - Vi-i), Vi = Vi{ui+i - Ui), i = 1, . . . ,m/2 (10) 

by grouping the variables into odd Ui = y2i-i and even Vi = y2i parts. Equation (fTUIl is bi- 
Hamiltonian with the quadratic Poisson bracket and Hamiltonian like the Volterra lattice |19j : 

{ui,Vi}o = UiVi {fj,nj+i}o = -Ujtij+i. (11) 

The Hamiltonians and the first integrals can be written in the new variables: 

m/2 ^ m/2 

Hl{u) = ^Ui + Vi, ffo(u) = - ^log(ui) +log(t'i), (12) 
j=l i=l 
m/2 ^ m/2 ^ 

■^9 = X! + ^i^i' Ic=^-ul + UiVi{Ui+Vi + Ui+i). (13) 

i=l ^ i=l 

Symplectic and time-reversible partitioned Runge-Kutta methods like the Lobatto IIIA-B 
methods |5j can be easily applied to Volterra lattice equations in the partitioned form. The 
symplectic Euler method which consists of a combination of explicit and implicit Euler is a first 
order Lobatto IIIA-B method. For the partitioned Volterra equations (jlUj) the symplectic Euler 
method becomes 

= + hu^iv^' - <+i) = + hv1+\u^_,, - <). (14) 

The second order Lobatto IIIA-B method for the partitioned Volterra lattice (|lflj) gives 

H = < (v^ + ^il - {vti + ^ili)) 1} = K + ^iDiv^+i - <) 

= + ^{kl + )) [v- + ^Ij - (Ci + ^^^i)) = K + (-^1 + \kli - K + )) 
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The internal stage vectors k"^ and are computed by solving a system of linear equations, 
whereas the vectors and P are obtained explicitly. The second order Lobatto IIIA-B method 
is known as Stormer-Verlet method for separable Hamiltonian systems. 

In order to show the preservation of the Poisson structure with the quadratic brackets 
we consider the corresponding two-form formulation [lijj 

duj Advj ^ duj A dvj+i 
ho ^^^^ ho ^^+1^^ 



Theorem 1 The Poisson structure of the Volterra lattice with the quadratic brackets (11) is 
preserved by the symplectic Euler method. 

Proof: In order to show the equality 

ho <^'-r' ^ho ~ho « ' 

we differentiate (11) to get 

dn^+i = dnr + /i(<+^ + (16) 
dv^+^ = d7;r + /i«+i-0'i<+^+K+^d<+i-H"^^^^<. (17) 

Using the equations ((TB|) and (fT7)) successively we obtain 

A = (1 + - <+i^))d< A - huyv'l+^ A 

d<+^Ad<+i = {l + hiv'l^,'-vr'))dvr'Adu^^,+hu?^,dv^^'Adv'l:;^,\ 

The equations (14) can be written in the equivalent form: 

,n+l 



1 + h{vr' - vt^) = 1 - h{ui^, - <) 

Combining all these into the two-forms in (jlSf) we obtain 

„,n+l n+l 

A = ' ^\ [dul A dv^ + hv^+^du'l A (iu^+i) - /in^d^.Y A , 



n+l n+l 

d<+iAd<+/ = ^±p^(d^^"Ad<+i-K+id<Adnr+i)+/i<+id<+iAd<+Y- 

Because du^Adu^j^^ = 0, the /i-order terms in the parenthesis vanish and taking in the summation 
over i and considering periodicity of the Volterra lattice one can easily see that the second terms 
of order h cancel and the quadratic Poisson brackets (11) are preserved by the symplectic Euler 
method. 

Because the second order Lobatto IIIA-B method is a composition of symplectic Euler meth- 
ods with step sizes /i/2, it preserves the quadratic Poisson bracket of the Volterra lattice. But 
higher order Lobatto IIIA-B methods can not preserve the Poisson structure, because they can 
not be written as combination of symplectic Euler method. Only diagonally implicit partitioned 
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Runge-Kutta methods can be written as combination of symplectic Euler method (see pp. 180, 

Recently finite dimensional systems with a Poisson structure arising after semi-discretization 
of certain partial differential equations have been integrated by splitting methods. An example 
of this is the preservation of the Lie-Poisson structure of Landau-Lifschitz lattice in partitioned 
form using a staggered scheme which corresponds to the second order Lobatto IIIA-B [S]. It was 
shown in |3j that the finite dimensional system which arises using a variational approximation of 
the time-dependent Schrodinger equation by Gaussian wave packets inherits a Poisson structure 
and various splitting methods were considered for its integration. 

The linear integrals are preserved exactly by all Runge-Kutta methods. The implicit mid- 
point rule preserves the quadratic integrals exactly assuming that the underlying system of 
nonlinear equations is solved within the machine accuracy. The Lobatto IIIA-B methods preserve 
only the quadratic integrals of the form Q = p^Dq, where D is an arbitrary matrix of appropriate 
dimension (Ch. 4, ,8J. Unfortunately the quadratic integral Iq (13) of the splitted Volterra 
lattice is not in this form. Higher order polynomial integrals like the cubic integral Ic and non- 
polynomial integrals like the Casimir function Hq are not preserved exactly by the symplectic 
Euler and Lobatto IIIA-B methods. 

For the periodic Volterra lattice we have used the following initial condition 

y{xi) = 1 + -^sech^ixi), Xj = -1 + (i - 1)-^-, i = l,...,m. 

All computations are done with a constant time step At = 0.1 over the time interval t G [0, 2000] 
for a Volterra lattice of dimension m = 20, 40, 80. The errors in the Hamiltonians and conserved 

quantities are given in Table 1 in the mean square root norm \Jj2iLi{I^ — I^Y ■, where P 
denote the computed Hamiltonians or first integrals at time step ti and N is the number of time 
steps. 

For all values of m and At, the Hamiltonian Hi is preserved with almost the same high 
accuracy for both methods. The Casimir Hq, the quadratic and cubic first integrals can not be 
preserved exactly, but the errors do not grow with time as in non-symplectic methods. Similar 
numerical results are obtained for the Toda lattice (see pp. 385-386 8^ and pp. 430-431 [2]). It 
was shown in |H1, Theorem 3.1, pp. 353, that for completely integrable systems, the symplectic 
integrators preserve the first integrals over long-time with an error 0{lSff)^ where p denotes the 
order of the method. The Lobatto IIIA-B method results in smaller errors than the symplectic 
Euler method because it is a second order accurate method. One can also observe that the 
Casimir Hq is preserved slightly better than the quadratic and cubic first integrals 
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Figure 1: Errors in the Hamiltonians and first integrals: symplectic Euler method, m = 40 
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Figure 2: Errors in the Hamiltonians and first integrals: Lobatto IIIA-B method, m = 40 
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Table 1: Average errors of the Hamiltonians and first integrals 







O V 111 Uit;^j Li^ J_JLlid llidjilUU. 


T nhatfn TTTA-R mpfhnr] 


m 


At 


Hi Hq Iq Ic 


Hi Hq Iq Ic 


20 
90 
20 


0.2 
n 1 
0.05 


2.589 -16 7.021 -08 6.200 -08 2.380 -07 
9 740 -16 1 940 -08 1 079 -08 4 100 -08 
1.652 -16 2.229 -09 1.853 -09 7.036 -09 


3.042 -16 2.698 -11 1.879 -10 1.597 -09 
9^96-16 5 157-19 968 -11 9 765 -10 
1.609 -16 9.373 -13 5.793 -12 4.896 -11 


40 
40 
40 


0.2 
0.1 
0.05 


3.016 -16 2.608 -10 2.158 -10 1.746 -09 
1.283 -16 4.820 -11 4.021 -11 3.218 -10 
2.309 -16 9.300 -12 7.863 -12 6.163 -11 


1.752 -16 7.557 -13 6.870 -12 6.036 -11 
4.824 -16 1.401 -13 1.202 -12 1.053 -11 
2.562 -16 2.540 -14 2.115 -13 1.857 -12 


80 
80 
80 


0.2 
0.1 
0.05 


2.903 -16 7.256 -12 7.162 -12 5.650 -11 
5.237 -16 1.314 -12 1.296 -12 1.017 -11 
6.117 -16 2.439 -13 2.396 -13 1.864 -12 


9.837 -16 1.409 -14 1.396 -13 1.229 -12 
2.244 -16 2.968 -15 2.342 -14 2.137 -13 
1.547 -16 6.175 -16 3.830 -15 3.745 -14 



4 Conclusion 

We have shown that the symplectic Euler method preserves the quadratic Poisson structure of 
the periodic Volterra lattice. The numerical results show excellent long time preservation of the 
Hamiltonian, Casimirs and the first integrals. Because of the singularity of the structure matrices 
we can obtain only local results for the backward error analysis in contrast to the global results 
obtained for Poisson systems with invertible structure matrices like in (Sec. IX. 3. 3, pp.297 jHj). 
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